final
Marcel Suleiman, Frank Pennekamp, Owen Petchey
7/27/2021
Introduction
This code contains all analyses done for the manuscript “Contrasting resistance and resilience of the coupled oxic and anoxic components of an experimental microbial community”. The first part deals with the oxygen data, the second part deals with the microbial 16S rRNA amplicons (PacBio full length sequencing)
Short summary of the project:
-8 micro-ecosystems, incubated at 24 °C for 8 days under light:dark cycles of 16:8h -afterwards, 4 micro-ecosystems (stressor treatment) were incubated in darkness for 7 days -afterwards, the stressor-treated columns were incubated again as the controls -Oxygen was measured at the top and bottom of the liquid part continuously every 5 min -samples were taken at day 8 (prior-stressor treatment), day 15 ( stressor-treatment), day 19 (short time recovery) and day 35 (long term recovery) at the height of the top and bottom oxygen sensor, respectively, for 16S rRNA full length amplicon sequencing
Aim of the study: Analysis of resistance and resilience of the oxygen concentration and the oxic/anoxic microbial community
Oxygen
data handling
oxygen_data <- read_excel("oxygen_data")Fix the date and times
oxygen_data <- oxygen_data %>%
mutate(Date_time = mdy_hms(paste0(Date, Time))) %>%
select(-Date, -Time)Calculate hourly average
hourly_value <- oxygen_data %>%
mutate(Date = date(Date_time),
Hour = hour(Date_time)) %>%
group_by(Sensor_Name, Date, Hour) %>%
dplyr::summarise(mean_Value = mean(Value)) %>%
mutate(Date_time = ymd_h(paste0(Date, Hour)))## `summarise()` has grouped output by 'Sensor_Name', 'Date'. You can override using the `.groups` argument.hourly_value <- hourly_value %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")))Calculate the DAILY mean, minimum, and maximum of the hourly means
daily_values <- hourly_value %>%
ungroup() %>%
group_by(Date, Sensor_Name) %>%
dplyr::summarise(mean_O2 = mean(mean_Value),
min_O2 = min(mean_Value),
max_O2 = max(mean_Value),
amplitude_O2 = max_O2 - min_O2) %>%
mutate(Date_time_midpoint = ymd_hms(paste0(Date, "12:00:00"))) %>%
pivot_longer(names_to = "Variable", values_to = "Oxygen", cols = 3:6)## `summarise()` has grouped output by 'Date'. You can override using the `.groups` argument.put in missing values for the days where there is no data
missing_data <- crossing(Sensor_Name = unique(daily_values$Sensor_Name),
Variable = unique(daily_values$Variable),
Date = dmy(paste(c(24:27), "12-2020")),
Date_time_midpoint = ymd_hms(paste0(Date, " 12:00:00")))%>%
mutate(Oxygen = NA)
daily_values <- rbind(daily_values, missing_data) %>%
mutate(Sensor_Name=tolower(Sensor_Name)) %>%
separate(Sensor_Name, into= c("Treatment", "Sensorposition","Column"),
remove=FALSE)
daily_values <- daily_values %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))plot the data
hourly values
plot_oxygen_a <- hourly_value %>%
mutate(Day2 = Day*24 + Hour) %>%
ggplot(aes(x = Day2, y = mean_Value, col = Sensor_Name)) +
geom_line() +
geom_rect(xmin = 195, xmax = 360,ymin = 0, ymax = 100, col = "lightblue", fill= "lightblue") +
geom_line() +
geom_rect(xmin = 510, xmax = 640,ymin = 0, ymax = 100, col = "white", fill= "white") +
scale_color_manual(values =c("Control_top_1" = "black", "Control_top_2" = "black", "Control_top_3" = "black", "Control_top_4" = "black","Disturbance_top_1" ="red","Disturbance_top_2" ="red","Disturbance_top_3" ="red","Disturbance_top_4" ="red","Control_bottom_1" = "grey", "Control_bottom_2" = "grey", "Control_bottom_3" = "grey", "Control_bottom_4" = "grey","Disturbance_bottom_1" ="grey","disturbance_bottom_2" ="grey","disturbance_bottom_3" ="grey","Disturbance_bottom_4" ="grey")) +
theme_bw() + theme(panel.grid = element_blank()) + scale_x_continuous(breaks = seq(0, 817, by=48), labels=seq(0,35, by=2))+
theme(legend.position = "none") +
ylab("Hourly mean oxygen value [%]") +
xlab("Day") +
geom_point(data=data.frame(Day2=c(190, 360, 456, 840), mean_Value=rep(-1, 4)), colour="black", size=3, shape=4)
plot(plot_oxygen_a)
just top sensor with all 4 variables and daily mean
plot_oxygen_b <- daily_values %>%
filter(Sensor_Name %in% c("control_top_1","control_top_2","control_top_3","control_top_4","disturbance_top_1","disturbance_top_2","disturbance_top_3","disturbance_top_4")) %>%
ggplot(aes(x = Day, y = Oxygen, col = Sensor_Name)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_line() +
facet_wrap( ~ Variable2, ncol = 4) +
scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), Oxygen=rep(-7, 4)), shape=4 ,colour="black", size=2) +
ylab("Daily mean oxygen value [%]") +
xlab("Day")
plot(plot_oxygen_b)
Applying t-test for every day to get p.value for top sensors
analysis <- daily_values %>%
filter(Sensor_Name != "control_top_1") %>%
na.omit() %>%
group_by(Date, Variable, Date_time_midpoint, Sensorposition) %>%
do(m1 = tidy(lm(Oxygen ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest() %>%
filter(term=="Treatmentdisturbance", Sensorposition == "top")## Warning: `cols` is now required when using unnest().
## Please use `cols = c(m1)`analysis <- analysis %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))Plot the estimates
plot_estimate_oxygen<- analysis %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_ribbon(aes(ymin=conf.low,ymax=conf.high, x=Day), alpha=.4)+
geom_rect(xmin = 22, xmax = 25,
ymin = -Inf, ymax = 100, col = "white", fill="white", alpha=1) +
facet_wrap( ~ Variable2, ncol = 4) +
scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), estimate=rep(-25, 4)), colour="black", size=2, shape=4) +
ylab("Estimate") +
xlab("Day")
print(plot_estimate_oxygen)
Microbial community
Load phyloseq data, remove chloroplast sequences
ps <- readRDS("ps_oxic_anoxic.rds")
ps_new <- subset_taxa(ps, Order != "Chloroplast")
taxa_names(ps_new) <- paste0("Seq", seq(ntaxa(ps)))Calculation of alpha diversity index
diversity_test<- estimate_richness(ps_new, split=TRUE, measures=NULL)## Warning in estimate_richness(ps_new, split = TRUE, measures = NULL): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.diversity_test$ID <- rownames(diversity_test)
names <- read_excel("barcodes_FKIII.xlsx")## New names:
## * `` -> ...1names1 <- as.data.frame(names)
names1$ID <- gsub("--", "..", names1$...1)
diversity_analysis <- merge(diversity_test, names1, by=c("ID"))
stderr <- function(x, na.rm=FALSE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
div_means <- diversity_analysis %>%
dplyr::filter(Position %in% c("bottom", "top")) %>%
dplyr::group_by(Position, Day, Treatment) %>%
dplyr::summarize(Shannon_mean = mean(Shannon),
Shannon_SE = stderr(Shannon))## `summarise()` has grouped output by 'Position', 'Day'. You can override using the `.groups` argument.div_raw <- diversity_analysis
plot_a_diversity_top <- ggplot() +
geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
geom_point(data = subset(div_raw, Position == "top"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
scale_shape_manual(values=c(17,16))+
geom_point(data=subset(div_means, Position == "top"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
geom_line(data=subset(div_means, Position =="top"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
ggtitle("top layer community") +
theme(plot.title = element_text(size = 11))## Warning: Ignoring unknown parameters: widthprint(plot_a_diversity_top)## Warning: Width not defined. Set with `position_dodge(width = ?)`
plot_a_diversity_bottom <- ggplot() +
geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
geom_point(data = subset(div_raw, Position == "bottom"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
scale_shape_manual(values=c(17,16))+
geom_point(data=subset(div_means, Position == "bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
geom_line(data=subset(div_means, Position =="bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
ggtitle("bottom layer community")+
theme(plot.title = element_text(size = 11))## Warning: Ignoring unknown parameters: widthprint(plot_a_diversity_bottom)## Warning: Width not defined. Set with `position_dodge(width = ?)`
estimates of the alpha diversity analysis
analysis_shannon <- div_raw %>%
group_by(Day, Position) %>%
do(m1 = tidy(lm(Shannon ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
# plot the estimates
plot_estimates_shannon_top <- analysis_shannon%>%
filter(Position=="top") %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day") +
scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35))
print(plot_estimates_shannon_top)
plot_estimates_shannon_bottom <- analysis_shannon%>%
filter(Position=="bottom") %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day") +
scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35))
print(plot_estimates_shannon_bottom)
#Anova
## top community
diversity_analysis_top<- diversity_analysis %>%
filter(Position=="top")
lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_top)
autoplot(lm_treatment)## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
anova(lm_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 0.0138 0.0138 0.0423 0.8385271
## Day 1 6.5483 6.5483 20.0874 0.0001225 ***
## Treatment:Day 1 1.6860 1.6860 5.1719 0.0311192 *
## Residuals 27 8.8018 0.3260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_top)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5282 -0.3754 0.1258 0.4857 0.6453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.51565 0.31188 17.685 2.25e-16 ***
## TreatmentDisturbance -0.84603 0.44627 -1.896 0.0687 .
## Day -0.06889 0.01441 -4.782 5.47e-05 ***
## TreatmentDisturbance:Day 0.04647 0.02044 2.274 0.0311 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.571 on 27 degrees of freedom
## Multiple R-squared: 0.4838, Adjusted R-squared: 0.4264
## F-statistic: 8.434 on 3 and 27 DF, p-value: 0.0004077##top sensor shannon index is affected "postponed", and it does not show resilience
Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_top)
autoplot(Anova_treatment)
anova(Anova_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 0.0138 0.0138 0.0988 0.756165
## as.factor(Day) 3 11.7919 3.9306 28.1259 7.133e-08 ***
## Treatment:as.factor(Day) 3 2.0299 0.6766 4.8418 0.009353 **
## Residuals 23 3.2143 0.1398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(Anova_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_top)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68415 -0.19311 0.01744 0.30636 0.44183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3096 0.1869 28.406 < 2e-16 ***
## TreatmentDisturbance -0.4739 0.2643 -1.793 0.08616 .
## as.factor(Day)15 -0.6177 0.2643 -2.337 0.02852 *
## as.factor(Day)19 -1.9468 0.2643 -7.365 1.72e-07 ***
## as.factor(Day)35 -1.9154 0.2643 -7.246 2.25e-07 ***
## TreatmentDisturbance:as.factor(Day)15 0.0423 0.3891 0.109 0.91436
## TreatmentDisturbance:as.factor(Day)19 0.8408 0.3738 2.249 0.03438 *
## TreatmentDisturbance:as.factor(Day)35 1.1886 0.3738 3.179 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3738 on 23 degrees of freedom
## Multiple R-squared: 0.8115, Adjusted R-squared: 0.7541
## F-statistic: 14.14 on 7 and 23 DF, p-value: 5.529e-07AIC(lm_treatment, Anova_treatment)## df AIC
## lm_treatment 5 58.94428
## Anova_treatment 9 35.71636#bottom community
diversity_analysis_bottom<- diversity_analysis %>%
filter(Position=="bottom")
##bottom sensor shannon index is affected directly, and it does show resilience
lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_bottom)
autoplot(lm_treatment)
anova(lm_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 6.0279 6.0279 7.9775 0.008631 **
## Day 1 22.0637 22.0637 29.1999 9.209e-06 ***
## Treatment:Day 1 0.2105 0.2105 0.2786 0.601757
## Residuals 28 21.1570 0.7556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_bottom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7771 -0.6599 0.0257 0.7631 1.4147
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.85624 0.47482 8.121 7.67e-09 ***
## TreatmentDisturbance 1.18319 0.67150 1.762 0.08898 .
## Day -0.07561 0.02193 -3.448 0.00181 **
## TreatmentDisturbance:Day -0.01637 0.03102 -0.528 0.60176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8693 on 28 degrees of freedom
## Multiple R-squared: 0.5722, Adjusted R-squared: 0.5264
## F-statistic: 12.49 on 3 and 28 DF, p-value: 2.307e-05#bottom sensor shannon index is affected directly, and it does show resilience
Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_bottom)
autoplot(Anova_treatment)
anova(Anova_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 6.0279 6.0279 12.050 0.001978 **
## as.factor(Day) 3 25.2738 8.4246 16.841 4.204e-06 ***
## Treatment:as.factor(Day) 3 6.1515 2.0505 4.099 0.017537 *
## Residuals 24 12.0059 0.5002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(Anova_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_bottom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22400 -0.31305 0.00521 0.39579 1.33089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.39590 0.35364 9.603 1.08e-09 ***
## TreatmentDisturbance -0.06541 0.50012 -0.131 0.897026
## as.factor(Day)15 -0.73628 0.50012 -1.472 0.153962
## as.factor(Day)19 -1.14224 0.50012 -2.284 0.031520 *
## as.factor(Day)35 -2.10233 0.50012 -4.204 0.000315 ***
## TreatmentDisturbance:as.factor(Day)15 1.83580 0.70728 2.596 0.015862 *
## TreatmentDisturbance:as.factor(Day)19 1.78258 0.70728 2.520 0.018783 *
## TreatmentDisturbance:as.factor(Day)35 0.11542 0.70728 0.163 0.871739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7073 on 24 degrees of freedom
## Multiple R-squared: 0.7573, Adjusted R-squared: 0.6865
## F-statistic: 10.7 on 7 and 24 DF, p-value: 4.674e-06AIC(lm_treatment, Anova_treatment)## df AIC
## lm_treatment 5 87.57164
## Anova_treatment 9 77.44136Calculation rel abundance, just sequences that appear > 0.1 %, making dataframes for each column position
ps_rel <- transform_sample_counts(ps_new, function(x) x / sum(x) )
relab_threshold <- 0.001
ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new)## [1] 9028ntaxa(ps_relab)## [1] 1856ps_relative <- transform_sample_counts(ps_relab, function(x) x / sum(x))
df_rel <- psmelt(ps_relative)
df_rel_top <- df_rel %>%
filter(Position == "top")
df_rel_bottom <- df_rel %>%
filter(Position == "bottom")NMDS
NMDS Plot
mds_whole <- ps_relative@otu_table %>%
as.data.frame() %>%
metaMDS(.,
distance = "bray",
k = 2, ## number of dimensions to reduce to
try = 200, ## number of random starts to try
autotransform = FALSE ## best not to use
)## Run 0 stress 0.1530756
## Run 1 stress 0.1863542
## Run 2 stress 0.1530827
## ... Procrustes: rmse 0.001494103 max resid 0.01029464
## Run 3 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 1.70895e-05 max resid 0.0001091102
## ... Similar to previous best
## Run 4 stress 0.2410315
## Run 5 stress 0.1571977
## Run 6 stress 0.1835709
## Run 7 stress 0.1571814
## Run 8 stress 0.2427357
## Run 9 stress 0.1530759
## ... Procrustes: rmse 0.0001919743 max resid 0.001269544
## ... Similar to previous best
## Run 10 stress 0.1530827
## ... Procrustes: rmse 0.001487175 max resid 0.01029609
## Run 11 stress 0.1530756
## ... Procrustes: rmse 1.930118e-05 max resid 0.0001257336
## ... Similar to previous best
## Run 12 stress 0.1530827
## ... Procrustes: rmse 0.001492498 max resid 0.01029447
## Run 13 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 6.204864e-06 max resid 3.648567e-05
## ... Similar to previous best
## Run 14 stress 0.1530827
## ... Procrustes: rmse 0.001486057 max resid 0.0102877
## Run 15 stress 0.1530829
## ... Procrustes: rmse 0.001510395 max resid 0.01031492
## Run 16 stress 0.221017
## Run 17 stress 0.1571977
## Run 18 stress 0.2474408
## Run 19 stress 0.1863542
## Run 20 stress 0.1864998
## Run 21 stress 0.1530829
## ... Procrustes: rmse 0.001516108 max resid 0.01031984
## Run 22 stress 0.1530829
## ... Procrustes: rmse 0.001514972 max resid 0.01031947
## Run 23 stress 0.1571977
## Run 24 stress 0.1530828
## ... Procrustes: rmse 0.001500227 max resid 0.01030554
## Run 25 stress 0.1530831
## ... Procrustes: rmse 0.001536167 max resid 0.01033908
## Run 26 stress 0.1530827
## ... Procrustes: rmse 0.001486437 max resid 0.01029374
## Run 27 stress 0.1571977
## Run 28 stress 0.153083
## ... Procrustes: rmse 0.001527379 max resid 0.01034268
## Run 29 stress 0.1530827
## ... Procrustes: rmse 0.001485112 max resid 0.01028405
## Run 30 stress 0.1530829
## ... Procrustes: rmse 0.001518242 max resid 0.01032105
## Run 31 stress 0.2354868
## Run 32 stress 0.1530827
## ... Procrustes: rmse 0.00148408 max resid 0.01027525
## Run 33 stress 0.1530756
## ... Procrustes: rmse 1.728505e-05 max resid 0.0001132889
## ... Similar to previous best
## Run 34 stress 0.2064158
## Run 35 stress 0.153083
## ... Procrustes: rmse 0.001529418 max resid 0.01033644
## Run 36 stress 0.1530756
## ... Procrustes: rmse 7.619846e-05 max resid 0.0005015251
## ... Similar to previous best
## Run 37 stress 0.1530827
## ... Procrustes: rmse 0.00148686 max resid 0.01029488
## Run 38 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 4.640423e-06 max resid 3.031545e-05
## ... Similar to previous best
## Run 39 stress 0.153083
## ... Procrustes: rmse 0.001520156 max resid 0.01032214
## Run 40 stress 0.2078723
## Run 41 stress 0.1530828
## ... Procrustes: rmse 0.001493384 max resid 0.01029707
## Run 42 stress 0.1530758
## ... Procrustes: rmse 0.0001460359 max resid 0.0009646334
## ... Similar to previous best
## Run 43 stress 0.1530827
## ... Procrustes: rmse 0.001485414 max resid 0.01028575
## Run 44 stress 0.1530827
## ... Procrustes: rmse 0.001483245 max resid 0.01027089
## Run 45 stress 0.2440261
## Run 46 stress 0.1571977
## Run 47 stress 0.2378509
## Run 48 stress 0.1530756
## ... Procrustes: rmse 9.933842e-06 max resid 6.319209e-05
## ... Similar to previous best
## Run 49 stress 0.153083
## ... Procrustes: rmse 0.001527391 max resid 0.01033342
## Run 50 stress 0.1571977
## Run 51 stress 0.1530832
## ... Procrustes: rmse 0.001539839 max resid 0.0103376
## Run 52 stress 0.1530756
## ... Procrustes: rmse 1.91514e-05 max resid 0.000125398
## ... Similar to previous best
## Run 53 stress 0.1571977
## Run 54 stress 0.1530827
## ... Procrustes: rmse 0.001487554 max resid 0.01028301
## Run 55 stress 0.2411633
## Run 56 stress 0.1530827
## ... Procrustes: rmse 0.001485677 max resid 0.0102878
## Run 57 stress 0.1530829
## ... Procrustes: rmse 0.001515751 max resid 0.01031523
## Run 58 stress 0.1530827
## ... Procrustes: rmse 0.001483622 max resid 0.01027334
## Run 59 stress 0.1530827
## ... Procrustes: rmse 0.001484592 max resid 0.01027958
## Run 60 stress 0.2064158
## Run 61 stress 0.2064158
## Run 62 stress 0.1530831
## ... Procrustes: rmse 0.001534772 max resid 0.010333
## Run 63 stress 0.1530827
## ... Procrustes: rmse 0.001481509 max resid 0.01026136
## Run 64 stress 0.1530827
## ... Procrustes: rmse 0.001483763 max resid 0.01027145
## Run 65 stress 0.1530827
## ... Procrustes: rmse 0.00148102 max resid 0.01025742
## Run 66 stress 0.1530756
## ... Procrustes: rmse 1.256547e-05 max resid 8.161904e-05
## ... Similar to previous best
## Run 67 stress 0.1571814
## Run 68 stress 0.1571977
## Run 69 stress 0.1530827
## ... Procrustes: rmse 0.001483074 max resid 0.01026654
## Run 70 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 4.041303e-06 max resid 1.748873e-05
## ... Similar to previous best
## Run 71 stress 0.153083
## ... Procrustes: rmse 0.001525471 max resid 0.01032985
## Run 72 stress 0.1571977
## Run 73 stress 0.1571814
## Run 74 stress 0.153083
## ... Procrustes: rmse 0.001527291 max resid 0.01033307
## Run 75 stress 0.153083
## ... Procrustes: rmse 0.001524731 max resid 0.01034017
## Run 76 stress 0.1530827
## ... Procrustes: rmse 0.001485807 max resid 0.0102793
## Run 77 stress 0.1571814
## Run 78 stress 0.1571977
## Run 79 stress 0.1835697
## Run 80 stress 0.2064158
## Run 81 stress 0.153083
## ... Procrustes: rmse 0.00152606 max resid 0.01033003
## Run 82 stress 0.1571977
## Run 83 stress 0.1530828
## ... Procrustes: rmse 0.0014955 max resid 0.010302
## Run 84 stress 0.1530827
## ... Procrustes: rmse 0.001487352 max resid 0.01027989
## Run 85 stress 0.1530829
## ... Procrustes: rmse 0.001512958 max resid 0.01031999
## Run 86 stress 0.1530756
## ... Procrustes: rmse 1.27091e-05 max resid 8.291818e-05
## ... Similar to previous best
## Run 87 stress 0.1571977
## Run 88 stress 0.1571977
## Run 89 stress 0.1530827
## ... Procrustes: rmse 0.001483973 max resid 0.01027375
## Run 90 stress 0.1530756
## ... Procrustes: rmse 1.26172e-05 max resid 7.745678e-05
## ... Similar to previous best
## Run 91 stress 0.153083
## ... Procrustes: rmse 0.001522445 max resid 0.01032893
## Run 92 stress 0.1530827
## ... Procrustes: rmse 0.001487952 max resid 0.01028765
## Run 93 stress 0.1530756
## ... Procrustes: rmse 8.25915e-06 max resid 5.378966e-05
## ... Similar to previous best
## Run 94 stress 0.1530827
## ... Procrustes: rmse 0.00148237 max resid 0.01026688
## Run 95 stress 0.1530827
## ... Procrustes: rmse 0.00148412 max resid 0.01027836
## Run 96 stress 0.2400825
## Run 97 stress 0.1530829
## ... Procrustes: rmse 0.001520601 max resid 0.01032439
## Run 98 stress 0.1571814
## Run 99 stress 0.153083
## ... Procrustes: rmse 0.001523792 max resid 0.010326
## Run 100 stress 0.2810977
## Run 101 stress 0.1571977
## Run 102 stress 0.1530828
## ... Procrustes: rmse 0.001507353 max resid 0.0103211
## Run 103 stress 0.1571977
## Run 104 stress 0.1571977
## Run 105 stress 0.1530829
## ... Procrustes: rmse 0.001510038 max resid 0.01031664
## Run 106 stress 0.1530827
## ... Procrustes: rmse 0.001484665 max resid 0.01028164
## Run 107 stress 0.1530756
## ... Procrustes: rmse 8.341181e-05 max resid 0.0005500282
## ... Similar to previous best
## Run 108 stress 0.1530829
## ... Procrustes: rmse 0.001512642 max resid 0.01031862
## Run 109 stress 0.2159959
## Run 110 stress 0.1530756
## ... Procrustes: rmse 9.08985e-06 max resid 5.798618e-05
## ... Similar to previous best
## Run 111 stress 0.1530756
## ... Procrustes: rmse 2.916728e-06 max resid 1.136089e-05
## ... Similar to previous best
## Run 112 stress 0.2323035
## Run 113 stress 0.1530831
## ... Procrustes: rmse 0.00154653 max resid 0.01035967
## Run 114 stress 0.1530827
## ... Procrustes: rmse 0.001481561 max resid 0.01025542
## Run 115 stress 0.1571814
## Run 116 stress 0.1530829
## ... Procrustes: rmse 0.001508675 max resid 0.01031415
## Run 117 stress 0.1571977
## Run 118 stress 0.1863542
## Run 119 stress 0.1530833
## ... Procrustes: rmse 0.001558055 max resid 0.01035742
## Run 120 stress 0.1530828
## ... Procrustes: rmse 0.001495171 max resid 0.01030318
## Run 121 stress 0.1530756
## ... Procrustes: rmse 9.641594e-06 max resid 6.317314e-05
## ... Similar to previous best
## Run 122 stress 0.1571977
## Run 123 stress 0.1530829
## ... Procrustes: rmse 0.001518184 max resid 0.01032365
## Run 124 stress 0.1530756
## ... Procrustes: rmse 6.272833e-05 max resid 0.0004042269
## ... Similar to previous best
## Run 125 stress 0.1863542
## Run 126 stress 0.1530756
## ... Procrustes: rmse 8.021191e-05 max resid 0.0005291918
## ... Similar to previous best
## Run 127 stress 0.2423328
## Run 128 stress 0.1530828
## ... Procrustes: rmse 0.001503851 max resid 0.01031026
## Run 129 stress 0.2363129
## Run 130 stress 0.1530828
## ... Procrustes: rmse 0.001505441 max resid 0.0103117
## Run 131 stress 0.1571977
## Run 132 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 3.263649e-06 max resid 1.890926e-05
## ... Similar to previous best
## Run 133 stress 0.1530827
## ... Procrustes: rmse 0.001483634 max resid 0.01027232
## Run 134 stress 0.1530756
## ... Procrustes: rmse 5.700027e-06 max resid 3.479409e-05
## ... Similar to previous best
## Run 135 stress 0.1530828
## ... Procrustes: rmse 0.001497463 max resid 0.01030401
## Run 136 stress 0.1530756
## ... Procrustes: rmse 1.136082e-05 max resid 7.484264e-05
## ... Similar to previous best
## Run 137 stress 0.1530827
## ... Procrustes: rmse 0.001482386 max resid 0.01026545
## Run 138 stress 0.2428061
## Run 139 stress 0.1530828
## ... Procrustes: rmse 0.001505543 max resid 0.01031146
## Run 140 stress 0.1530756
## ... Procrustes: rmse 5.186104e-05 max resid 0.0003416616
## ... Similar to previous best
## Run 141 stress 0.1571814
## Run 142 stress 0.1530827
## ... Procrustes: rmse 0.00148759 max resid 0.01029645
## Run 143 stress 0.1530827
## ... Procrustes: rmse 0.001483407 max resid 0.01027115
## Run 144 stress 0.1571977
## Run 145 stress 0.1835717
## Run 146 stress 0.2215934
## Run 147 stress 0.1530756
## ... Procrustes: rmse 4.645247e-06 max resid 2.494342e-05
## ... Similar to previous best
## Run 148 stress 0.1530756
## ... Procrustes: rmse 8.509436e-06 max resid 5.423649e-05
## ... Similar to previous best
## Run 149 stress 0.1530756
## ... Procrustes: rmse 5.541056e-06 max resid 3.638999e-05
## ... Similar to previous best
## Run 150 stress 0.1571814
## Run 151 stress 0.1530827
## ... Procrustes: rmse 0.001488455 max resid 0.01029962
## Run 152 stress 0.1530829
## ... Procrustes: rmse 0.001515513 max resid 0.01032507
## Run 153 stress 0.1530829
## ... Procrustes: rmse 0.001517771 max resid 0.01032287
## Run 154 stress 0.153083
## ... Procrustes: rmse 0.001527008 max resid 0.01032672
## Run 155 stress 0.1530756
## ... Procrustes: rmse 1.001727e-05 max resid 6.597879e-05
## ... Similar to previous best
## Run 156 stress 0.1530831
## ... Procrustes: rmse 0.001538281 max resid 0.01034477
## Run 157 stress 0.1530756
## ... Procrustes: rmse 4.027731e-05 max resid 0.0002651449
## ... Similar to previous best
## Run 158 stress 0.1530756
## ... Procrustes: rmse 6.453689e-06 max resid 4.171174e-05
## ... Similar to previous best
## Run 159 stress 0.1571814
## Run 160 stress 0.153083
## ... Procrustes: rmse 0.001521279 max resid 0.01032701
## Run 161 stress 0.2356746
## Run 162 stress 0.1530827
## ... Procrustes: rmse 0.001489317 max resid 0.0103119
## Run 163 stress 0.1530829
## ... Procrustes: rmse 0.001512308 max resid 0.01030889
## Run 164 stress 0.153083
## ... Procrustes: rmse 0.001526984 max resid 0.01033047
## Run 165 stress 0.1571977
## Run 166 stress 0.1571814
## Run 167 stress 0.1530828
## ... Procrustes: rmse 0.001498739 max resid 0.01030449
## Run 168 stress 0.1530828
## ... Procrustes: rmse 0.001496279 max resid 0.01030118
## Run 169 stress 0.2044065
## Run 170 stress 0.2395851
## Run 171 stress 0.1530756
## ... Procrustes: rmse 4.440006e-06 max resid 2.539756e-05
## ... Similar to previous best
## Run 172 stress 0.2044065
## Run 173 stress 0.1530829
## ... Procrustes: rmse 0.001514033 max resid 0.01032076
## Run 174 stress 0.1530828
## ... Procrustes: rmse 0.001495232 max resid 0.0102853
## Run 175 stress 0.1530756
## ... Procrustes: rmse 7.904381e-06 max resid 5.194586e-05
## ... Similar to previous best
## Run 176 stress 0.1530831
## ... Procrustes: rmse 0.001546085 max resid 0.01035948
## Run 177 stress 0.1530829
## ... Procrustes: rmse 0.001514981 max resid 0.01031598
## Run 178 stress 0.2078904
## Run 179 stress 0.1530831
## ... Procrustes: rmse 0.001539807 max resid 0.01034065
## Run 180 stress 0.1530831
## ... Procrustes: rmse 0.001549926 max resid 0.01036804
## Run 181 stress 0.1571977
## Run 182 stress 0.1530827
## ... Procrustes: rmse 0.001482756 max resid 0.01026846
## Run 183 stress 0.1530758
## ... Procrustes: rmse 0.0001462399 max resid 0.0009590674
## ... Similar to previous best
## Run 184 stress 0.1530828
## ... Procrustes: rmse 0.001499389 max resid 0.01030584
## Run 185 stress 0.1530827
## ... Procrustes: rmse 0.001484218 max resid 0.01027437
## Run 186 stress 0.1530827
## ... Procrustes: rmse 0.001485497 max resid 0.01028503
## Run 187 stress 0.1530827
## ... Procrustes: rmse 0.001482626 max resid 0.01026788
## Run 188 stress 0.1530828
## ... Procrustes: rmse 0.00149758 max resid 0.01030293
## Run 189 stress 0.1530756
## ... Procrustes: rmse 6.892817e-06 max resid 3.665205e-05
## ... Similar to previous best
## Run 190 stress 0.1530829
## ... Procrustes: rmse 0.001521049 max resid 0.01032749
## Run 191 stress 0.4002155
## Run 192 stress 0.1530827
## ... Procrustes: rmse 0.001488543 max resid 0.0102806
## Run 193 stress 0.1530827
## ... Procrustes: rmse 0.001485343 max resid 0.01028032
## Run 194 stress 0.1530829
## ... Procrustes: rmse 0.001510011 max resid 0.01031812
## Run 195 stress 0.2429624
## Run 196 stress 0.1530847
## ... Procrustes: rmse 0.001636921 max resid 0.01049326
## Run 197 stress 0.1571977
## Run 198 stress 0.1571977
## Run 199 stress 0.1530828
## ... Procrustes: rmse 0.001498565 max resid 0.01030433
## Run 200 stress 0.1530756
## ... Procrustes: rmse 2.143592e-05 max resid 0.0001398277
## ... Similar to previous best
## *** Solution reached## 0.16
mds_whole_res <- ps_relative @sam_data %>%
as.tibble() %>%
select(Treatment, Column, Day, Position,Day_name) %>%
bind_cols(as.tibble(scores(mds_whole, display = "sites")))## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectnmds_plot_supplements <- ggplot(mds_whole_res, aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(shape = Position, color=Treatment),
size = 2) +
scale_shape_manual(values =c(17,16))+
facet_wrap("Day", ncol=2)+
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))
print(nmds_plot_supplements)
NMDS separated for NMDS1 and NMDS2 axis
NMDS_top <- mds_whole_res %>%
filter(Position=="top")
NMDS_bottom <- mds_whole_res %>%
filter(Position=="bottom")
# Dynamics
plot_NMDS1_top <- NMDS_top %>%
ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=16, size = 2)+
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("top layer community") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS1_top)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS2_top <- NMDS_top %>%
ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=16, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
ggtitle("top layer community")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS2_top)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS1_bottom <- NMDS_bottom %>%
ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=17, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("bottom layer community")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS1_bottom)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS2_bottom<- NMDS_bottom %>%
ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=17, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("bottom layer community")+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS2_bottom)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
estimate analysis for NMDS data
analysis_NMDS1_top <- NMDS_top %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS1_top <- analysis_NMDS1_top %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS1_top)
analysis_NMDS2_top <- NMDS_top %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS2_top <- analysis_NMDS2_top %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS2_top)
analysis_NMDS1_bottom <- NMDS_bottom %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS1_bottom <- analysis_NMDS1_bottom%>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS1_bottom)
analysis_NMDS2_bottom <- NMDS_bottom %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS2_bottom <- analysis_NMDS2_bottom%>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS2_bottom)
final plot of figures combined
p1 <- plot_oxygen_a
p2 <- plot_oxygen_b
p3 <- plot_a_diversity_top
p4<- plot_a_diversity_bottom
p5 <- plot_estimates_shannon_top
p6 <- plot_estimates_shannon_bottom
p7 <- plot_NMDS1_bottom
p8 <- plot_NMDS2_bottom
p9 <- plot_NMDS1_top
p10 <- plot_NMDS2_top
p11 <- plot_estimate_oxygen
p12<- plot_estimates_shannon_top
p13<- plot_estimates_shannon_bottom
p14 <- plot_estimate_NMDS1_bottom
p15 <- plot_estimate_NMDS2_bottom
p16 <- plot_estimate_NMDS1_top
p17 <- plot_estimate_NMDS2_top
#Figure 1
plot_1 <- p1/
p2/
p11
plot_1 +
plot_layout(heights=c(2,1,1)) + plot_annotation(tag_levels = "a")
#Figure 2
plot <- (p4|p13)/
(p3|p12)/
(p7|p14)/
(p8|p15)/
(p9 |p16)/
(p10 |p17)
plot + plot_annotation(tag_levels = "a")## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Width not defined. Set with `position_dodge(width = ?)`## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
#Supplementary Information
Visualization of rel.abundance for supplementary information
df_class<- df_rel %>%
group_by(Class, Treatment,Day, Column, Position) %>%
dplyr::summarize(Abundance=sum(Abundance))## `summarise()` has grouped output by 'Class', 'Treatment', 'Day', 'Column'. You can override using the `.groups` argument.index <- which(df_class$Abundance>=0.05)
class_to_keep <- unique(df_class[index,"Class"])
class_to_keep <- unname(unlist(class_to_keep))
df_class$Class_filter <- ifelse(df_class$Class %in% class_to_keep, df_class$Class,"Zother")
cols <-c("Actinobacteria"="bisque1", "Alphaproteobacteria" = "yellow", "Babeliae"="lightcyan1","Bacilli" = "green4", "Bacteroidia" ="maroon1", "Campylobacteria" ="royalblue1","Zother"="black","Chlorobia" ="thistle1", "Clostridia"= "firebrick1", "Cyanobacteriia" = "magenta4", "Fimbriimonadia" = "grey","Gammaproteobacteria"="red", "Microgenomatia"="grey38","Oligoflexia"="snow1", "Phycisphaerae"="cornflowerblue", "Planctomycetes"="lightskyblue1","Spirochaetia"="lavender", "Vampirivibronia" = "lightsalmon", "Verrucomicrobiae" = "slateblue3")
supplement_plot1 <-df_class%>%
ggplot(aes_string(x = "Column", y = "Abundance", fill="Class_filter" )) +
geom_bar(stat = "identity", position = "stack", col="black") +
facet_wrap(Position~Day, ncol = 4) +
theme(axis.text.x = element_text(angle = -45, hjust = 0, vjust=0.5, size=9))+
ylab("Relative abundance") +
scale_fill_manual("legend", values = cols)
print(supplement_plot1)